#include "AnaDef.h"
#include "style.h"

#define EPSILON 1e-10

void drawMCTOFpurity(const Int_t tofid, const Bool_t ktpconly)
{
  style::SetGlobalStyle(0);

  TFile *_file0 = TFile::Open("/u/xlu/.task/Fragmentation/TrainOutput/TOF13_MCjet_23022014/MC10de.root");

  const Int_t tpccutmode = 1;
  const TString dn=Form("PWGJE_FragmentationFunction_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
  const TString ln=Form("fracfunc_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);    

  TDirectory *fdir=_file0->GetDirectory(dn);
  if(!fdir){
    printf("no fdir!!\n");exit(1);
  }

  TList * ll=(TList*)fdir->Get(ln);
  if(!ll){
    printf("no list!!\n");exit(1);
  }

  THnSparseD* hh=(THnSparseD*)ll->FindObject("THnIDFF");
  for(int ii=0; ii<hh->GetNdimensions(); ii++){
    printf("%d: %s\n", ii, hh->GetAxis(ii)->GetTitle());
  }
  /*
1: 3
0: TrackEta
1: JetPt
2: TrackTPCsig
3: Log10TrackP
4: Log10TrackPt
5: z
6: xi
7: pdg
8: comb
9: tof
10: tpc
   */

  TAxis * axj=hh->GetAxis(1);
  axj->SetRangeUser(10.001, 15-0.001);

  TAxis * axpt=hh->GetAxis(4);
  axpt->SetRangeUser(TMath::Log10(0.3), TMath::Log10(3));


  TH2D *h2=hh->Projection(ktpconly?10:9,7);
  h2->GetXaxis()->SetRangeUser(0,3);
  h2->GetYaxis()->SetRangeUser(0,3);
  style::ResetStyle(h2);

  TCanvas *c1=new TCanvas;
  style::PadSetup(c1);
  c1->SetRightMargin(0.03);
  gStyle->SetOptTitle(1);

  gStyle->SetOptStat("");
  //
  Double_t inclusion=-999;
  Double_t exclusion=-999;
  if(!ktpconly){
    if(tofid==1){
      inclusion = 3;
      exclusion = 2.5;
    }
    else if(tofid==2){
      inclusion = 3;
      exclusion = 3;
    }
    else if(tofid==3){
      inclusion = 3;
      exclusion = 3.5;
    }
    else{
      printf("bad tofid %d\n", tofid); exit(1);
    }
  }
  //

  //renormaliz rows
  for(Int_t iy=1; iy<=h2->GetNbinsY(); iy++){
    const Double_t nor = h2->Integral(h2->GetXaxis()->FindBin(0.), h2->GetXaxis()->FindBin(3.), iy, iy);
    if(nor<EPSILON)
      continue;

    for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
      const Double_t ent = h2->GetBinContent(ix, iy);
      h2->SetBinContent(ix,iy, ent/nor);
      h2->SetBinError(ix,iy, TMath::Sqrt(ent)/nor);
    }
  }

  h2->SetTitle(Form("0: p, 1: #pi, 2: K, 3: e (jet pT(GeV/c) %.0f - %.0f %s) %.2e < p_{T} GeV/c < %.2e;MC true type;%s", axj->GetBinLowEdge(axj->GetFirst()), axj->GetBinUpEdge(axj->GetLast()) , ktpconly?"TPC 5][5]":Form("TOF %.1f][%.1f", inclusion, exclusion), pow(10,axpt->GetBinLowEdge(axpt->GetFirst())), pow(10,axpt->GetBinUpEdge(axpt->GetLast())), ktpconly?"PreID":"TOF"));
  h2->SetMarkerSize(3);
  h2->SetLineWidth(2);
  h2->Scale(100);
  h2->SetMinimum(1);
  h2->Draw("box text");
  c1->Print(Form("MCtofpurityTOF%d_ktpconly%d.eps", tofid, ktpconly));
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=3){
    printf("wrong argc!!\n"); return 1;
  }

  drawMCTOFpurity(atoi(argv[1]), atoi(argv[2]));

  return 0;
}


